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Abstract — The natural trajectory tracking problem is 
studied for generic quantum states represented by density 
operators. A control design based on the Hilbert-Schmidt 
distance as a Lyapunov function is considered. The control 
dynamics is redefined on an extended space where the 
LaSalle invariance principle can be correctly applied 
even for non-stationary target states. LaSalle's invariance 
principle is used to derive a general characterization of 
the invariant set, which is shown to always contain the 
critical points of the Lyapunov function. Critical point 
analysis of the latter is used to show that, for generic 
states, it is a Morse function with n\ isolated critical points, 
including one global minimum, one global maximum and 
n\ — 2 saddles. It is also shown, however, that the actual 
dynamics of the system is not a gradient flow, and therefore 
a full eigenvalue analysis of the linearized dynamics about 
the critical points of the dynamical system is necessary to 
ascertain stability of the critical points. This analysis shows 
that a generic target state is locally asymptotically stable 
if the linearized system is controllable and the invariant 
set is regular, and in fact convergence to the target state 
(trajectory) in this case is almost global in that the stable 
manifolds of all other critical points form a subset of 
measure zero of the state space. On the other hand, if 
either of these sufficient conditions is not satisfied, the 
target state ceases to be asymptotically stable, a center 
manifold emerges around the target state, and the control 
design ceases to be effective. 

Index Terms — trajectory tracking, LaSalle invariance 
principle, eigenvalue analysis, quantum systems. 

I. Introduction 

Recent advances in quantum optics, trapping of 
cold atoms, ions and molecules, and breakthroughs 
in nano-scale engineering of artificial atoms and 
molecules have prompted significant interest in 
ways to control these systems and the development 
of the theoretical foundations of quantum control 
theory. One of the major concerns is how to design 
the dynamics of a given system to steer its state to 



a desired target and stabilize it in this state. One 
proposed technique is Lyapunov control, where the 
control function is designed such that a suitably cho- 
sen Lyapunov function is monotonically decreasing 
along every trajectory. A number of Lyapunov- 
function-based control designs have been proposed 
and numerous results established (see, e.g., [ 1 ] — [ 11] 
and references therein). The convergence analysis in 
most of these works is based on the application of 
LaSalle's invariance principle. However, the invari- 
ant set for quantum systems is usually large and thus 
convergence to the target state cannot be inferred 
from the invariance principle directly. 

For pure-state quantum systems with states de- 
scribed by wave-functions, the setting considered 
in most of the papers to date, there are many 
results. For example, it was shown in [6] that a 
particular variant of the method is effective under 
certain sufficient conditions, equivalent to control- 
lability of the linearized system. [7] also proposed 
a modified control design based on an "implicit" 
Lyapunov function for systems that do not satisfy 
the conditions for local asymptotic stability. Al- 
though pure states play a crucial role in quantum 
mechanics, they form a set of measure zero at the 
boundary of the domain of density operators, and 
in practice quantum systems are often not in pure 
states to begin with, due to imperfect preparation, 
or because the system is an ensemble of many 
quantum particles such as atoms or molecules. For 
this reason it is important to consider control in 
the context of density operators representing generic 
quantum states. It is also important to consider non- 
stationary target states, i.e., target states that are not 
eigenstates of the system Hamiltonian, because it 
is these so-called superposition states that exhibit 
non-classical behavior such as interference and en- 
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tanglement, which is a crucial resource for novel 
quantum applications such as quantum information 
processing. 

This more general case was considered more 
recently in [10], where the convergence properties 
of the control design were investigated for both pure 
and mixed target states, including non- stationary 
states. Simulations for bilinear control systems sug- 
gest that under certain strict sufficient conditions as 
given in [10], all trajectories starting outside the 
invariant set converge to the target state as opposed 
to other points in the invariant set. Simulations also 
suggest that when these conditions are not satisfied 
then we converge to states (trajectories) in the 
vicinity of the target, yet never reach the target, and 
the states (trajectories) we appear to converge to in 
the latter case are not critical points of the Lyapunov 
function. It has been suggested that this could be 
due to numerical errors, but this fails to explain 
why such errors should not affect the convergence 
equally in all cases. Can we find analytical results 
that explain these observations? Another issue is 
the correct application the invariance principle for 
non- stationary target states. Simply transforming to 
a rotating frame to make the target state formally 
appear stationary, as has been done in several pa- 
pers, is problematic because the resulting system 
is non- autonomous, making the application of the 
invariance principle problematic. 

To address these issues and evaluate possible 
explanations e.g., that the invariant set consists only 
of critical points of the dynamical system, which 
are repulsive except for the target state, necessitates 
several steps, including careful characterization of 
the invariant set, the set of critical points of the 
Lyapunov function, and a stability analysis of the 
critical points of the dynamical system. Analysis 
of the critical points of the Lyapunov function as 
considered, e.g., in [18] recently, is not sufficient, 
as the dynamics is not a gradient flow. However, we 
can still use this analysis to show that the Lyapunov 
function is a Morse function if we restrict our 
attention to generic states. This allows us to draw 
certain conclusions about the dimensions of the 
stable (unstable) manifolds of the hyperbolic critical 
points of the dynamical system, from which we can 
draw general conclusions about the effectiveness of 
the method in steering the system to the target. 

The paper is organized as follows. In Sec. [XT] 
the control problem is defined and some basic 



issues such as different notions of convergence 
for non- stationary target state are briefly discussed. 



In Sec. Ill we formulate the control dynamics as 
an autonomous dynamical system defined on an 
extended space including the system and the target 
state, allowing us to apply LaSalle's invariance prin- 
ciple [15] to obtain a characterization of the LaSalle 
invariant set. This set is shown to depend on the 
Hamiltonian as well as the target state. In Sec. IV 
we determine the critical points of the Lyapunov 
function on the extended space and analyze their sta- 
bility for generic target states. In Sec. [V] we analyze 
the stability of the critical points of the dynamical 
system in terms of the eigenvalues of the linearized 
system. Rigorous results are derived for generic 
stationary target states and somewhat weaker sta- 
bility results are given in the non- stationary case. 



In Section VI it is shown explicitly that relaxing 
the strict requirements on the Hamiltonian leads to 
loss of asymptotic stability due to the emergence 
of center manifolds about the target state. Finally, 
we compare and analyze the difference between our 
result and the argument in one preceding work [10]. 

II. Natural Trajectory Tracking Problem 
for Quantum Systems 

A. Quantum states and evolution 

The state of a quantum system defined on an n- 
dimensional Hilbert space H ~ C n can be repre- 
sented by a density operator p, i.e., antixn positive 
hermitian operator with unit trace, whose evolution 
is governed by the following equation: 

ihp(t) = [H,p(t))=Hp(t)-p(t)H, 

where H is an nxn Hermitian operator representing 
the system Hamiltonian, and we shall choose units 
such that h = 1. When the system is not closed, i.e., 
interacts with an external environment, additional 
terms are required to account for dissipative effects, 
and in the Markovian case the evolution is described 
by Lindblad master equation [13], for example. 
In the following, we will restrict our analysis to 
Hamiltonian systems and to an important class of 
mixed states, we shall refer to as generic states. The 
same analysis can be applied to density operators 
representing non-generic states, although the results 
will be different. 

Definition 11.1. A density operator p represents a 
generic state if it has n distinct eigenvalues. 



3 



B. Control problem 

We study the bilinear Hamiltonian control system 

p(t) = -i[H + f(t)H 1 ,p(t)}, (1) 

where f(t) is an admissible real-valued control field 
and H and H 1 are a system and control interaction 
Hamiltonian, respectively, both of which will be as- 
sumed to be time-independent. The general control 
problem is to design a certain control function f(t) 
such that the system state p{t) with p(0) = po 
converges to the target state p d for t — > oo. We 
shall assume here that the initial and target states, 
po and p d , have the same spectrum because this 
is a necessary condition for the target state to be 
reachable under unitary evolution. 

Since the free Hamiltonian H can usually not be 
turned off, any target state will evolve according to 

Pd (t) = -i[H , p d (t)}. (2) 

It is easy to see that p d is stationary if and only if 
it commutes with H , [H , p d (0)} = 0. For any state 
that does not commute with H the quantum control 
problem becomes a "natural tracking problem" [14], 
where the objective generally is to find a control 
f(t) such that the trajectory p(t) with initial state 
po under the controlled evolution asymptotically 
converges to the trajectory of p d (t). 

C. Trajectory tracking vs orbit tracking 

When the target state p{t) is non-stationary, we 
can have two different control objectives. We could 
require p{t) — > p d (t), as t — > oo, which is known as 
trajectory tracking; alternatively, we could require 
p(t) — > 0(p d (t)), which is known orbit tracking, 
where 0(p d (t)) is the orbit of p d , defined to be the 
set of points p d (t) passes through. By definition, 
the notion of orbit tracking is weaker than that of 
trajectory tracking. In particular, the set 0(p d (t)) 
can be very large for states that follow non-periodic 
trajectories, which is the case for most states, except 
for systems of Hilbert space dimension 2, and higher 
dimensional systems for which the eigenvalues of 
Hq are commensurate (rational multiples of each 
other). Therefore in this work we focus on quantum 
state control in the sense of trajectory tracking as 
this is the strongest notion of convergence and well- 
defined for any trajectory. 



III. Lyapunov-based Control 

A natural design for the control f(t) is in- 
spired by the construction of a Lyapunov func- 
tion V(p,p d ). We try to find a control such that 
^ t V(p(t), p d (t)) < 0, i.e., V decreases along any 
trajectory. If V(p(t), p d (t)) decreases to zero, we 
have pit) -> p d {t). 

Define M. to be the set of density operators 
isospectral with p d (0). M. is a compact manifold, 
whose dimension depends on the spectrum of p d (0). 
For a generic p d with n distinct eigenvalues we 
have M ~ U(n)/U(l) x ... x U(l), where the 
denominator has n factors and M. has dimension 
n 2 — n. Consider the joint dynamics for (p(t), p d (t)) 
on .M x M: 

p(t) = -i[H Q + f(t)H l ,p(t)}, (3a) 
p d (t) = -i[H , p d (t)}. (3b) 

The Hilbert-Schmidt norm \\A\\ = ^Tr(A^A) in- 
duces a natural distance function on M. x M., 
which provides a natural candidate for a Lyapunov 
function 

V(p, Pd ) = \\\P-P4 2 = ^Tr[(p-p d ) 2 ]. (4) 

Since p and p d are required to be isospectral, we 
have Tr(p 2 ) = Tr(p^), and hence 

V(p,p d )=Tr[p 2 d (t)]-TT[p(t)p d (t)], (5) 

the Lyapunov function used in [10]. If p d = \ijj d ){'4' d \ 
and p = \ip){ip\ we have furthermore 

v(ipM = i-\{Mt)\m)\ 2 , (6) 

a Lyapunov function often used for pure-state con- 
trol [7]. Choosing 

f(p,p d ) = KTr([-iH 1 ,p]p d ), k>0, (7) 

we find V(p(t), p d (t)) < 0. Without loss of gener- 
ality, we set k — 1 in the following. 

Given the target state p d (0), the dynamics under 
the Lyapunov control is described by the following 
nonlinear autonomous dynamical system on M x 
M: 

p(t) = -t[H + f(p,p d )H u p(t)], (8a) 

pi(t) = -i[H ,p d (t)], (8b) 
f(p,p d )=Tr([-iH 1 ,p]p d ), (8c) 
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A. LaSalle invariance principle 

Theorem IIL1 (LaSalle invariance principle [15]). 

x = f(x) be an autonomous dynamical system, 
V(x) be a Lyapunov function on the phase space 
Q = {x}, satisfying V(x) > for all x ^ x 
and V(x) < 0, and let 0(x(t)) be the orbit of 
x(t) in the phase space. Then the invariant set 
E = {(D(x(t))\V(x(t)) = 0} contains the positive 
limiting sets of all bounded solutions, i.e., any 
bounded solution converges to E as t — > +oo. 

Remark III.l. The theorem holds for both real and 
complex dynamical systems, and what has basically 
been proved is that all bounded solutions with 
V(x) 7^ converge to the set of solutions with 
V(x) = 0. It does not matter if V(x) = for many 
x, as is the case for ([§]), for which V vanishes on 
the entire set {(p, pd) G Ai x Ai : p = pd}. 

The quantum system ([8]) defined on the extended 
phase space Ai x Ai is autonomous and any solution 
(p(t), pd(t)) is bounded. Although the Lyapunov 
function ([5]) is not positive definite, V — if and 
only if p = pd, which is sufficient to apply the 



LaSalle invariance principle III.l to obtain: 



Theorem III.2. Any trajectory (p(t),pd{i)) under 
the Lyapunov control ([7]) will converge to the invari- 
ant set E = {(p 1; p 2 ) e M x M : V(p(t), p d (t)) = 
0, (p(0),p d (0)) = ( Pl ,p 2 )}. 

We note here that except when pd is a stationary 
state, we must consider the dynamical system on the 
extended phase space Ai x Ai since V(p, pd) is not 
well-defined on Ai . Having established convergence 
to the LaSalle invariant set E, the next step is to 
characterize E for the dynamical system ([8]). 

B. Characterization of the LaSalle Invariant Set 

LaSalle's invariance principle reduces the con- 
vergence analysis to calculating the invariant set. 
It has been argued (e.g. [10], [11]) that the in- 
variant set consists of all points p that commute 
with the target state pd, i.e., [pd, p] = 0. However, 
this characterization is only valid for ideal systems 
and stationary target states. Thus we shall first 
reconsider the characterization of the invariant set. 
Following standard techniques in nonlinear stability 
analysis [12], we shall see that the invariant set of 
the autonomous dynamical system ([8]) defined on 
the extended state space Ai x Ai comprises all 



pairs (pi,p2) whose commutator is diagonal even 
for ideal systems, and it is much larger for non- 
ideal systems. 

It is easy to see that V = is equivalent to f(t) = 
and a standard argument shows that 

= f = Tr([-iH 1 ,p]p d ) 

= f = Tr([-iH u p]p d ) + Tt{[-iH 1} p\p d ) 

= -Tr([Ad_ iHo (-i£fi), p]p d ) 







where B m 



Tr([Adl^ Htf^pK), 



Ad r 



a AHo (—iHi) represents m-fold 
commutator adjoint action of — iH on —iH\. Not- 
ing Tr([A,B]C) = Tr([A, [B, C\\) shows that the 
LaSalle invariant set consists of all (pi, p 2 ) such that 



Ti(B m [ Pl , p 2 ]) = 0. 
= span{S m }^=f and B S Q 



(9) 



oo 



If we set B s 

span{5 m }™=o w i tn Bq = —iU\ then ^ is equiv 
alent to [p, p d ] orthogonal to the subspace B = B\ 
with respect to the Hilbert-Schmidt inner prod- 
uct. The Lie algebra su(n) can be decomposed 
into an abelian part, the Cartan subalgebra C = 
span{Afc}fc=i, and an orthogonal subspace T, which 
is a direct sum of n(n — l)/2 root spaces spanned 
by pairs of generators {A/^, A^} (see Appendix 

Choosing a Hilbert space basis such that H is 
diagonal H = diag(a l5 . . . , a n ) with a k > a k+1 
and uju — ag — a k , which is always possible as H 
is Hermitian, setting bu = ceke + ifiki and expanding 
—iHi G su(n) with respect to the basis (25) gives 



71—1 



fe=i 



£=k+l 



Bo, 



2m- 1 



B 



2m 



Noting that H is diagonal and using (27) gives 

71—1 71 

E E( 

k=i e=k+i 

71—1 71 

E E (~^-) mu kT[ a k£^ke — Pk£^k£\ 
k=l t=k+l 

(10) 

Definition III.l. The Hamiltonian of the dynamical 
system ([§]) is called ideal if 

(i) H is strongly regular, i.e., ^ uj pq unless 
(k,£) = (p,q). 

(ii) Hi is fully connected, i.e., bki ^ except 
(possibly) for k — I. 
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Theorem III.3. The subspace B n ~ n generated by 
the Ad-brackets is equal to T if and only if the 
Hamiltonian is ideal. 

Proof: The technique in this proof is a direct 
application of the property of Vandermonde matrix, 
which has also been applied to discuss the con- 
trollability [16], [17]. Since the dimension of T is 
n 2 — n and B m £ T for all m > 0, it suffices to 
show that the elements B m for m = 1, . . . , n 2 — n 
are linearly independent. Moreover, the subspaces 
spanned by the odd and even order elements, 



B s 



odd 



span{i?2m-i : 1 < 2m — 1 < s} and 
span{I?2m : 1 < 2m < s}, respectively, 



are orthogonal since 

r> r> / i \m+m 

-D2m-1-D2m' — {— J- J 



2m— 1 2m' v 



— 3 f J(6 fc £)Q ; (6 fc /^)A fe £A fc /£/ + Q : (6^)3 f ?(6fc'F)A^Afe^'], 
and thus observing the equalities 

Tr(A^Afc/£/) = Tr(Afc£Afc/^/) = — 25 kk '5w (11a) 
Tr(A H A^) = (Hb) 
shows that for all m, m! > 



Tr( J B 2m _ 1 S 2m = (-1) 



m+m 



'EE- 



2m- 1, ,2m' 



UJ 



k'f 



ft(b M )Z(b ki )[-\ 2 ke + \ 2 H ] = 0. 



odd 



Thus it suffices to show that the elements of B\ 

2 

and i3™ ve ~™ are linearly independent separately. 

For the odd terms, suppose there exists a vector 
c = (ci, . . . , c s ) T of length s = n(n — l)/2 such 
that Em=i c mJ B 2m _i = 0. Noting that u kk = and 
{ujik) 2 = (—u k e) 2 this gives n{n — l)/2 non-trivial 
equations 

s 



-a; 



2 \m-l 



m=l 



(12) 

for 1 < k < £ < n. Since iQ k t ^ 0, &^ 7^ by 
hypothesis, Eq. (12) can be reduced to VLc = 0, 
where is a matrix: 



n 



-UJ 



12 



0J 



12 



-6J 



2 \m-l 
12 J 



13 



UJ 



13 



-6J 



2 \m-l 



13; 



-6U, 



\m— 1 



n— l,n/ 



3) 



Since is a Vandermonde matrix, condition (ii) 



of the proposition guarantees that Eq. (12) has 
only the trivial solution c = 0, thus establishing 
linear independence. For the even terms we obtain 
a similar system of equations, which completes the 
proof. ■ 

Theorem III.4. Assuming the Hamiltonian o/([S]) is 
ideal, (/?i,p 2 ) belongs to the invariant set E if and 
only if [pi, p 2 ] = diag(ci, . . . , c„). 

Proof: The necessary part follows from Eq. (|9]) 
and Theorem III. 3 By Q the states in the in- 
variant set have to satisfy Tr(B[p, pj) = for 
all B e B n -™, and by Theorem III.4 we have 
B n ~ n = T for ideal systems. Thus, [p, p<f] must be 
in the subspace orthogonal to T, which is the Cartan 
subspace, i.e., the diagonal (trace-zero) matrices. 
For the sufficient part note that p k (t) = 
1,2, 



e~ iHot p k e iHot , k 



-iH t 



pie 



iH t 



-iH t e iH t 



-iH t 



\PuP2\e 

and e~ lHot diagonal. Thus if [pi,p 2 ] 



iH t 



diag(ci, . . . ,c„ 
diag(ci, . . . ,c v 



then 



-iH t 



iH t 



[pi,p 2 ] and hence (pi,p 2 ) G E. 



Thus we have fully characterized the invariant 
set for systems with ideal Hamiltonians. The re- 
sult also shows that even under the most stringent 
assumptions on the Hamiltonian, the invariant set 
is generally much larger than the desired solution. 
Therefore, the invariance principle alone is not 
sufficient to establish convergence to the target state. 



Remark III.2. Theorem III. 3 shows that for a 
system with ideal Hamiltonian the Ad-brackets span 
the entire tangent space to the state manifold, 

2 

B n ~ n = T, and thus the linearization defined 
on the tangent space is controllable. It is also 
easy to verify that this condition is necessary for 
controllability of the linearization for generic states. 

IV. Lyapunov Function as Morse Function 

We show the LaSalle invariant set of ([8]) always 
contains the critical points of the Lyapunov function 

V, and characterize the stability of these critical 
points. Notice that the stability here refers to the 
stability of these critical points as stationary states 
of the gradient flow induced by the Lyapunov 
function. The gradient flow in general need not be 
related to the particular dynamics of the system, 
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e.g., prescribed by the equation of motion ([8]). 
Indeed we shall see that ([8]) is not a gradient 
flow. The stability analysis of the critical points of 
V is still useful, however, because under certain 
conditions, for instance, when the stationary state 
of ([8]) is hyperbolic, the dimensions of the stable 
(unstable) manifold of ([8]) agree with those of the 
stable (unstable) manifold of the gradient flow. 

We start with the case where pd is given, and 
V(p, pa) is effectively a function of p on M. Similar 
topics have also been discussed in [18], but the 
major point here is that we can find the critical 
points of V(pi,p 2 ) on M. x M. from the critical 
points of V(p) on M.. Since p can be written as 
p = U p d W for some U in the special unitary group 
§U(n), V can also be considered as a function on 
SU(n), V(U) = V(Up d U ] pd). It is easy to see that 
the critical points of V(p) correspond to those of 
V(U), and as Tr(p^) is a constant for a given p d , it 
is equivalent to find the critical points U E §U(n) 
of 

J{U) = Tr(p 2 d ) - V(U) = Tx{Up d tfp d ). (14) 
Lemma IV. 1. The critical points Uq of J(U) defined 



by (14) are such that [po, pd] = for p = UqpJJq. 



Proof: Let {<7 m }m=i be the orthonormal basis 
of su(n) given in Appendix |A| Any U E SU(ra) near 
the identity / can be written as U = e x ' a , where 
a = (<ti, . . . , cr n 2_ 1 ) and x E W l is the coordinate 
of U, and any U in the neighborhood of Uq can be 
parameterized as U = e x ' a UQ. Thus (14) becomes 



J = Tr[(e x - d U )p d (U f e- S - S )p, 



(15) 



At the critical point U , V J = implies 

= Tr{a m [U oPd Ul p d \) Vm. (16) 

Thus \UqpJUq, pd] E su(n) is orthogonal to all basis 
elements o m , and therefore [U oPd UQ, p d \ = 0. ■ 
More generally, for V( P x, p 2 ) defined on M. x A4, 
with pi = U\p d U\ and p 2 = U 2 p d U 2 , the critical 
points of J(pi,p 2 ) = Tr(p!p 2 ) and V(p 1 ,p 2 ) = 
Tr(p 2 ) - Tr(pxp 2 ) coincide as Tr(p^) = Tr(p^) is 
constant. 

J(Ux, U 2 ) = Tr (u xPd U\U 2 p d U\ 

= Tr {{UlU^pdiUlU^t 



together with Lemma |IV.1| shows that J attains its 
critical value when [{U 2 Ui)p d (U 2 Ui)\ p d ] = and 

= U 2 [(U 2 j U 1 )pd(U j 2 U i y,pd]U j 2 
= \U lP dUlU 2 pdUl} = \p u p 2 }. 

Thus we have the following: 

Theorem IV. 1. For a given target state p d (0), the 
critical points of the Lyapunov function V(pi,p 2 ) 
onMxM satisfy [pi, p 2 ] = and thus belong to 
the LaSalle invariant set E. 

Next, for a given generic p d , J(p) = Tr (ppd), and 
thus V(p) = V(p,pd), are Morse functions on AL 
i.e. its critical points are hyperbolic [21]: 

Theorem IV.2. If p d is generic then J(p) = Tr (ppd) 
is a Morse function on Ai. Two of the n! hyperbolic 
critical points correspond to the global maximum 
and minimum of J, respectively, and the other n\ — 2 
points are saddles with critical values Jo satisfying 

Jm'm ^ Jo ^ Jma,x- 

Proof: See [20] for the definition of sink, 
source and saddles. For a given p d , p is a critical 
point of J(p) if and only if [po, p d ] = 0, i.e. there 
exists a basis such that 



Pd 
Po 



diag(u>i, . . 
diag(uv ( i), 



where r is a permutation of the numbers {1, . . . , n}. 
Since pd is generic and hence the Wk are distinct, 
there are n\ distinct permutations and thus n! critical 
points with critical values J(po) = ELi^^ffc)' 

Next, in order to calculate the Hessian matrix, we 
need to find a parameterization of the points near p . 
Recalling that for any p E M. we have p = U pdU^ 
for some U E SU(n), consider J as a function on 
SU(n) with J(U) = Tr (Up d Wp d ). Let U be a 
critical point of J(U) and p = U p d UQ. Any U 
in the neighborhood of Uq can be parameterized 
as U = e x ' a Uo, where x E MJ 1 _1 and a is the 
orthonormal basis for su(n) defined in appendix [A] 
Substituting this into J gives 

J(x) = Ti(e^U p d U^e-^pd) 

= Tr[(I + x ■ a + \{x ■ a) 2 ) U p d U^ x 

(I-x -a+ \{x-a) 2 )p d ] + 6(|x| 3 ) 
= Tr[C/ p d Ul p d ] + Tr[(f ■ B) 2 popd] 
- Tr[(x • a)p (x ■ a)p d ] + 6(|£| 3 ), 
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where we have used p pd = PdPo- Taking the basis 
to be a =JX k , X ki , Xki} with X k , X ke and X kji as in 
we find that the Hessian matrix „ 9 i 



appendix A 

at p is diagonal, i.e., the basis vectors are eigenvec- 
tors. The first n — 1 diagonal entries corresponding 
to X k vanish but as we are only interested in the 
tangent space to the manifold spanned by {A^, Xke}, 
we can restrict our attention to this subspace. On 
this subspace, i.e., for tjj = Xu or X k e, we have 

d 2 J 
d 2 Xj 



2Ti[a 2 p p] - 2 Tr {a jpoCTjp}. 



The action of a* 



>i — Xke or Oj — am 
the subspace spanned by the basis vectors e k and 
e£. On this subspace A| £ is identity operator and 
the conjugate action of a a on the diagonal matrix 
p swaps its fc-th and £-th diagonal entries. Since 
p is non-degenerate, any swap X k e or X k e will 
make either larger or smaller than zero. Thus 
the Hessian matrix at p is diagonal with n 2 — n 
non-zero diagonal entries, corresponding to n 2 — n 
independent directions in the tangent space of M.. 
Therefore, all n\ critical points p are hyperbolic, 
and J is a Morse function. The maximal critical 
value occurs only when p = pd and the minimal 
value occurs only when w r (fc)' s are in an increasing 
order. For all other critical values, there always 
exists a swap that will increase the value of J 
and one that will decrease it, showing that they are 
saddle points. ■ 

V. Effectiveness of Lyapunov control for 

IDEAL SYSTEMS 

When the hyperbolic critical points of the Lya- 
punov function are also the stationary points of 
the dynamics, there are restrictions on the possible 
dynamics near those critical points. In particular, 
if the dynamics is the gradient flow of the the 
Lyapunov function then there is a simple correspon- 
dence between the number of negative (positive) 
eigenvalues at the critical point and the dimension 
of the stable (unstable) manifold at the critical point 
as a stationary solution. However, in general, this 
does not hold for a dynamical system other than the 
gradient flow. To be a gradient flow, the coefficient 
matrix of the linearized system has to be symmetric 
and we will see that ([8]) is not the gradient flow 
of any function. Therefore, in order to investigate 
the stability and to calculate the dimension of the 



X k p is restricted to 



stable manifold at any stationary point, we have to 
resort to the definition of the stable manifold, and 
investigate the linearized dynamics. 

Throughout this section we shall assume that the 
Hamiltonian is ideal. Without loss of generality we 
further assume that H has zero trace, as the identity 
part of H only changes the global phase. Once 
the Hamiltonian is chosen, the LaSalle invariant set 
E depends only on the target state p&. As stated 
before, throughout this paper we focus on generic 
states, i.e., assuming p d has n distinct eigenvalues, 
and assume p(0) and pd(0) have the same spectrum. 
Similar tools can be applied to non-generic states 
but they must be separately investigated as the 
topology of the critical points for non-generic states 
is different. 

A. Stationary (generic) target state 

We work in a basis where H is diagonal. If p d is 
stationary, i.e., [Ho, pd] = 0, then it is also diagonal, 
and ([8]) reduces to a dynamical system on M.: 

p(t) = -i[H + f{p)H 1 ,p(t)] (17a) 
f(p)=Tr{[-iH 1 ,p(t)]p d ) (17b) 

with the corresponding LaSalle invariant set 

E = {p :V(p(t)) = 0,p(0) = p } 

= {po ■ [po, Pd] = diag(ci, . . . , c n )} (18) 

Let pd = diag(u7i, . . . , w n ), with w k ^ We for k ^ L 
For any p e E, [pd, p] is diagonal if and only if p is 
diagonal, with diagonal elements as a permutation 
of (wi, . . . ,w n ). According to previous section, 
these n! stationary points are also the hyperbolic 
critical points of the Lyapunov function V(p). 

Theorem V.l. If pd is a generic stationary state 
then the invariant set contains exactly the n\ critical 

(k) 

points of the Lyapunov function p d , k = 1, . . . , n\, 
that commute with p d and have the same spectrum. 

These n\ points are the only stationary solutions 
and all the other solutions must converge to one 
of these points. Through analyzing the sign of 
the eigenvalues of the coefficient matrix of the 
linearized system, we shall see that p d is asymp- 
totically stable, and all other stationary points are 
unstable. In order to achieve this, we require a real 
representation for ( [T7j ). A natural choice is the Bloch 
representation. Let {^ k } k=1 be an orthonormal basis 
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for all n x n Hermitian matrices, with £„2 = ^1. 
We have p = EfcSfc£fc, with s k = Tr(p£ fc ). Since 
the dynamics is trace-preserving, i.e. s n 2 = ^= 
is constant, we can further reduce the dynamics 
onto the subspace 
as a vector s G 1 
dynamical system 



on — 1 



, and p can be represented 
nZ ~ 1 . Accordingly, the quantum 
([8]) can be represented as 



k{t) = {A Q + f(s,s d )A x )s{t) 
s d (t) = A s d (t) 
f(s, s d ) = s'/Axs, 

where A and A x are two anti-symmetric matrices: 



A (m,n) = Tr(iH Q [£ m ,£ n ]) 
Ai(m,n) = Tr(iHi[£ m ,£ n ]) 



When p d is stationary, this system reduces to 



k(t) = (A + f^A^t) 
f(s) = s d T Ais, 



(20) 
(21) 

(22a) 
(22b) 



and the Lyapunov function ([4]) is repres ented as 
V(s) = hWs — SdW 2 . According to Theorem |V. 1 1 for 



— Srf|| 2 . According to Theorem 
a generic p d , (22) has n\ stationary points, denoted 
as s^ k \ k = 1, . . . ,n\. The linearized system near 
the stationary state ^ is 



(23) 



where D f (^) = A +sfA 1 s<^A 1 +A 1 s< k ^.s d T A 1 



is a linear map defined on 



s/A^ 



0, and D f (^ 



i n -\ f(^) = gives 
= A + A 1 ^-s d T A 1 . 



Since A and A\ are anti- symmetric, Dj(s^ >) can- 
not be a symmetric matrix, and the dynamics cannot 
be a gradient flow of any function. Therefore, the 
topology near ^ as a critical point of V is not 
enough to infer the local dynamics in its vicinity, 
and we need to actually calculate the eigenvalues 
of D f (^). 



Remark V.l. The state space S_m of {22) is the set 



of all Bloch vectors s G IR n corresponding to 
density operators p G M.. For generic p d the state 
manifold M. is a flag manifold homeomorphic to 
§U(n)/ exp(C), where C is the Cartan subspace of 
the Lie algebra su(n) and exp(C) is its exponential 
image in §U(n), corresponding to diagonal unitary 
matrices with determinant 1. Hence, the tangent 
space Tm(po) of M. at any point po corresponds 
to the non-Cartan subspace T of su(n) and the 
Cartan elements of s\i{n) correspond to the tangent 



space of the isotropy subgroup of po- In the real 
representation, W 1 _1 15 therefore the direct sum of 
the (n 2 — n)-dimensional tangent space St to the 
manifold Sm and the (n— l)-dimensional subspace 
Sc corresponding to the Cartan subspace of su(n). 

Theorem V.2. For a generic stationary target state 
p d all the n\ stationary states of the dynamical 



system ( |77p are hyperbolic, i.e. all eigenvalues of 
Df(s^), restricted on St, have nonzero real parts, 
for k = 1, . . . , n\. Among those stationary states, p d 
is the only sink, all other points are saddles, except 
the global maximum, which is a source. 

Proof: Let s De one of the n\ stationary states. 
We first show that Df(so) vanishes on the [n — 1)- 
dimensional subspace Sc, which is orthogonal to 
St- In the second step we show that Df(s ) is 
invariant on St and has n 2 —n non-zero eigenvalues. 
Finally, we show that the restriction of Df(s ) onto 
St does not have any purely imaginary eigenval- 
ues, from which it follows that s is a hyperbolic 
stationary state, and the local dynamics of ( fT7| ) near 
every stationary state can therefore be approximated 
by the linearized system [20]. 

Lemma V.L Df(s ) vanishes on the subspace Sc- 

This lemma shows that sq is not a hyperbolic 
fixed point of the dynamical system ( [22] ) defined 
on M. n ~ 1 . However, we are only interested in the 
dynamics on the manifold Sm, an d thus it suffices 
to show that s is a hyperbolic fixed point of the 
restriction of Df(s ) to the tangent space St of S M . 

Lemma V.2. The restriction ofDf(s ) to St is well- 
defined, represented by a matrix B with n 2 — n non- 
zero eigenvalues. 

Lemma V.3. If i/3 is a purely imaginary eigenvalue 
of B then it must be an eigenvalue of B , i.e., 
ij3 = ±iuJke. for some (k, £), and either the asso- 
ciated eigenvector e must be an eigenvector of B 
with the same eigenvalue, or the restriction of A\s$ 
to the (k,£) subspace must vanish. 



Lemma V.3| shows that B can have a purely 
imaginary eigenvalue i/3 only if i(3 = ±iu>ke f° r 
some (k, £), and either u^*^ = 0, i.e., the projection 
of A\Sq onto the (k,£) subspace vanishes, or the 
associated eigenvector is also an eigenvector of B . 
In the first case this means that Ais vanishes on 
the subspace T ke , or equivalently that [—iHi, p ] has 
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no support in Tm, which contradicts the assumption 
that Hi is fully connected and p has non-degenerate 
eigenvalues. On the other hand, if e is an eigenvector 
of B with eigenvalue i(3 = ±itOke and H is 
strongly regular then the projection of e onto the 
(k,£) subspace is proportional to (l,±i) and e is 
zero elsewhere, and thus iFe = implies v^ k ' e ' = 0, 
which contradicts the fact that the projection AiSd or 
[—iHi,pd] onto the (k,£) subspace must not vanish 
if Hi is fully connected and pd has non-degenerate 
eigenvalues. Thus we can conclude that if H is 
strongly regular, Hi fully connected and p d has non- 
degenerate eigenvalues, Df(s ) cannot have purely 
imaginary eigenvalues, and thus s is hyperbolic. ■ 
This theorem illustrates that the n\ critical points 

P 



of V are also the 



nl 



states of (17). Since V(p a 



hyperbolic stationary 
= 0, p = p d must be 
a dynamical sink, with all eigenvalues of Df(sd) 
having negative real parts. Similarly, p = pW with 
V(p^) = Vmax must be a dynamical source. All 
the other p (fc) with < V^(p (fc) ) < V max must 
be saddles, with eigenvalues of Df(^) having 
both negative and positive real parts, for otherwise 
p( fc ) would be a sink or source, and thus a local 
minimum or maximum of V, which would contra- 
dict Theorem |IV.2[ Moreover, the dimension of the 
stable (unstable) manifold at p^ must agree with 
the the index number of V(p^) at p^ k \ This is 
a very useful observation as it allows us to infer 
that the dimension of the stable (unstable) manifold 
at p( fe ) is independent of the specific value of pd, 
only dependent on the relative location of p^ as a 
critical point of V and the system dimension n. The 
theorem also shows that each of the n\ — 2 saddle 
points, p( fe ) has a stable manifold. Solutions on the 
stable manifold converge to p^ and thus the saddles 
are not repulsive, as asserted in [10], and we can in 
construct counter-examples to Theorem 1 in [10]. 

For example, consider a three-level system with 
Hq strongly regular and Hi off-diagonal and fully 
connected. For a generic stationary target state such 
as p d = \ diag(3, 2, 1), the LaSalle invariant set 



6 stationary states — p 



(i) 



p(k) re f erre( i to as the antipodal 



Pd 



consists of 3! 
and five other 

points in [10]. The coefficient matrix Df(p^) of 
the linearized system has eigenvalues with negative 
real parts for every p( k > except the global maximum 



P 



(6) 



| diag(l, 2, 3), and thus four of the antipodal 
points have stable manifolds and solutions converg- 
ing to them. An even easier way to see that these 



points cannot all be repulsive is to note that if, e.g., 
p = \ diag(2, 3, 1) was repulsive then we would 
have V(p(t)) < V(p ) for all p(t) in a neighborhood 
of p , and thus p would be a local maximum of 
the Lyapunov V(p) = \ Tr(p — p d ) 2 , contradicting 
the fact that it is a saddle point of V. On the 
other hand, we note that any state p(t) starting 
outside the invariant set E has at least one off- 
diagonal component, and as Hi is fully connected 
and pd non-degenerate, the off-diagonal components 
of [— iHi, pd] are all nonzero. Thus, the trajectories 
converging to the saddle points satisfy conditions (1) 
and (2) of Theorem 1 in [10]. Condition (3) is also 
satisfied as CardjF t ([— iH\, pd]) = 3 = dim.M/2, 
where dim M. = 3 2 — 3. Thus by Theorem 1 in [10] 
they should converge to pd, which is not the case. 

Nonetheless, the stable manifolds of the unstable 
stationary states are not a serious obstruction to con- 
vergence. In fact, since all solutions not converging 
to pd are located on the union of the n\ — 2 stable 
manifolds of dimension < n 2 — n, which form a 
measure-zero set in the state space, we can conclude 
that almost all initial states converge to p d , i.e., p d 
can be considered almost globally asymptotically 
stable, and the Lyapunov design effective in this 
case. This is illustrated in Fig. [TJa) which shows 
that for a stationary generic pd, all simulated non- 
stationary trajectories with random p(0) converge to 
Pd exponentially. 

B. Generic non-stationary target state 

In this case characterizing the invariant set is 
more complicated as E may contain points with 
nonzero diagonal commutators. 

Example V.l. Let (p(0), p d (0)) = (pi,p2) with 



Pi 



j_ 

12 

_ J_ 
12 

12 



_ J_ 
12 
11 
24 
1 
8 



_J_ 
12 
1 
8 
11 
24 



P2 



12 



12 

1 
3 



12 



12 



Pi and P2 are isospectral and the commutator 
[pi,P2] = lli diag(0, 1, — 1) is diagonal, and thus 
(pi,p 2 ) G E, but [pi,p 2 ] ^ 0. 

When Pd(0) is chosen such that E contains 
points with nonzero diagonal commutators, Fig. |TJd) 
shows that all trajectories generated by the simu- 
lations fail to converge to pd(t), and the original 
control design becomes ineffective, even for systems 
with ideal Hamiltonians. Fortunately, however, the 
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Fig. 1: Time evolution of V(p(t), Pd{t)) with y-axis in logarithmic scale. Each graph shows V(p(t), Pd(t)) 
for N = 50 different initial states p(0). The graphs represent four different types of generic pd, of which 
(a,b,d) are for ideal Hamiltonian, and (c) for non-ideal Hamiltonian. (a) shows that for stationary pd, all 
trajectories converge exponentially to the target state to within machine precision. The negative slopes in 
(b) suggest that for a non- stationary target state with regular E, all simulated trajectories still converge 
to the target trajectory albeit at a slower rate compared to (a). For a non- stationary p d with irregular 
E as in (d), or a stationary p d with H\ not fully connected as in (c), on the other hand, the slopes of 
V(p(t), pd(t)) in the log-plot vanish at different finite distances from the target state for all simulated 
trajectories, indicating convergence to states or trajectories at various non-zero distances from the target. 




5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 



(a) Stationary target state (b) Non-stationary target state with regular E 




above example is quite exceptional. We shall see 
that E = {[pi, P2} = 0} still holds for a very large 
class of generic target states Pd(t), and in these cases 
Lyapunov control tends to be effective. For conve- 
nience, E is called regular if it only contains points 
with zero commutators, and irregular otherwise. 

Noting that we can write [pi, P2] = — Ad P2 (pi), 
where Ad P2 is a linear map from the Hermitian or 
anti-Hermitian matrices into su(n), let A(s 2 ) be the 
real [n 2 — 1) x [n 2 — 1) matrix corresponding to the 



Bloch representation of Ad P2 . Recall su(n) = T®C 
and W 2 - 1 = S r ® S c , where S c and S T are the 
real subspaces corresponding to the Cartan and non- 
Cartan subspaces, C and T, respectively. Let A(s 2 ) 
be the first n 2 — n rows of A(s 2 ) (whose image is 
SV). We have the following lemma, with proof in 
Appendix [Cj 

Lemma V.4. Given a generic pd(t), the LaSalle 
invariant set E is irregular if and only if 

rankyl(srf(0)) < n 2 — n. 
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This lemma provides a necessary and sufficient 
condition on Pd(0) to ensure that [pi, P2] diagonal 
implies [p 1; p 2 ] = 0. Assuming the first n 2 — n 
rows correspond to St, let Ai be the submatrix 
generated from the first n 2 — n rows and last 
n 2 - n columns of A(s d (0)). If det(Ai) ^ then 
rankA(s d (0)) = n 2 — n, hence E = {[pi,/^] = 0}. 
We can easily verify that if the_ diagonal elements of 
p d (0) are not equal then det(A 1 )det[A(p d (0))] = 0. 
is a non-trivial polynomial, i.e., det(Ai) can only 
have a finite set of zeros. Hence: 

Theorem V.3. The LaSalle invariant set for a 
generic pdif) is irregular only if det[A(pd(0))] = 
or some of the diagonal elements of Pd(0) are equal. 

Thus for most generic non-stationary p d (0) we 
still have E regular. In this case, given pd(0), 
let p d k \o) be the n\ critical points of V{p) = 
V(p,pd(0)) with critical values V^. Then we can 
easily see that the n! flows (p d k \t) , p d (t)) starting 
from (pj^(0), Prf(0)) with / = are solutions of the 
dynamical system satisfying [p d k \t) , pd(t)} = for 
any t, and thus the corresponding trajectories are in 
the LaSalle invariant set and are the critical points 
of V with V(p d k \t), p d (t)) = Vk, and we can show 
that any p(t) must converge to one of these critical 
trajectories. The trajectories with p d k \t) 7^ Pd{t) 
cannot be asymptotically stable as they correspond 
to unstable critical points of V. Furthermore, let 
Vk be the critical values of V ordered in an in- 
creasing sequence with V = 0, corresponding to 
the global minimum. Then all initial states p(0) 
with V(p(0), Pd(0)) < V\ must converge to p d (t) 
as V is monotonically decreasing, and thus pd(t) 
is locally asymptotically stable. We can summarize 
these findings in the following: 

Theorem V.4. Given a generic non-stationary pd(t) 
if the LaSalle invariant set E is regular, then any 
trajectory p(t) converges to one of the n\ critical 

1 n\. All critical tra- 



trajectories p d k \t), k = 1, 
jectories are unstable, except pdif), which is lo- 
cally asymptotically stable, and the global maximum 
p( n '\t), which is repulsive. 

As illustrated in Fig. [jjb), for regular E, all 
trajectories p(t) in computer simulations keep con- 
verging to p d (t). Although the convergence speed is 
slow compared to the case of generic stationary p d , 
the case of regular E is qualitatively different from 



the irregular one in Fig [TJd), where the rate of con- 
vergence drops to zero after some time, resulting in 
flat-lining of the trajectories in the semi-logarithmic 
plot. We conclude from these simulations that for 
a generic non- stationary pd(t) Lyapunov control is 
still effective even E is regular, although the conver- 
gence speed may be slow, while when E is irregular, 
the original control design becomes ineffective even 
for systems with ideal Hamiltonians. 

VI. (Non)Effectiveness of Lyapunov 
Control for Non-ideal Systems 

In the previous section we showed that p{t) = 
pd(t) is the only locally asymptotically stable tra- 
jectory if the system is ideal and the target state 
is regular. Realistic systems, unfortunately, often 
do not satisfy the strong Hamiltonian requirements, 
and we now show that in this case the target state, 
even if it is stationary, ceases to be a hyperbolic 
critical point. A center manifold emerges and most 
solutions do not converge to the target state, render- 
ing the method ineffective. This clearly illustrates 
that the dynamics ([8]) is very different from the 
gradient flow of the Lyapunov function and shows 
that the LaSalle invariance principle and critical 
point analysis of the Lyapunov function do not 
suffice to analyze the stability for realistic systems, 
and eigenvalue analysis of the linearized system is 
necessary. To fully understand the dynamics in this 
situation we need to analyze it case by case. In 
the following we present an analysis for three-level 
systems, which illustrates the techniques that can be 
applied to other cases. 

A. Hi not fully connected 

Assume H still strongly regular but H t be not 
fully connected, for example, consider 



Hn 



a 1 







a 2 







<*3 



Hi 





K 








0" 





where we assume a\ < a 2 < 03 and 61, 62 7^ 0. 
According to the characterization of the LaSalle 



invariant set E in Section III a necessary condi- 
tion for (pi, p 2 ) to be in the invariant set E is 
[pi,p 2 ] to be orthogonal to the subspace spanned 
by B = spanfAj^o with B m = Ad^-iHj. 
Comparison with (jTOb shows that if the coefficient 
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bk£ = then none of the generators B rn have 
support in the root space T ke of the Lie algebra, 
and it is easy to see that the subspace of su(n) 
generated by B is the direct sum of all root spaces 
7fc£ with bki 7^ 0. Thus, in our example, a necessary 
condition for (pi, p 2 ) to be in the invariant set E is 
[Pi-, P2] ^ ^13 ®C, which shows that we must have 



"ii a 13 
ct 22 
a* 3 a 33 



(24) 



Furthermore, if (pi, p 2 ) is of type (24) then 



W)[pi,p 2 ]W) 



"ii 




a: 



13 



«22 




e iuJl3t a 13 


«33 



with Uq = e~ lHot and u ki = at — a k , also has this 
form. Therefore, [pi,p 2 ] G C © T\ 3 is a necessary 
and sufficient condition for the invariant set E. For 
stationary generic pd, E consists of all (pi, P2) with 
p 2 = pa and 



Pi 



P11 (3 l3 
/3 22 




ft: 



13 



33 



Thus, the invariant set _E contains 3! = 6 stationary 
states corresponding to fii 3 = 0, which coincide 
with the critical points of V(p), and an infinite 
number of trajectories with /3 13 7^ 0. 

We check the stability of linearized system near 
these fixed points, concentrating on the local behav- 
ior near Sd. Working with a real representation of the 



linearized system (23) and using the same notation 
as before, we can still show that Df(sd) has n 2 — n 
nonzero eigenvalues, with n = 3 in our case. Since 
—iH\ has no support in the root space T\ 3 , the A13 
and A13 components of AiSd, (which correspond 
to [—iH^pd]) vanish, and Df(s^) has a pair of 
purely imaginary eigenvalues, whose eigenspaces 
span the root space T\ 3 , and four eigenvalues with 
non-zero real parts, which must be negative as Sd 
is locally stable from the Lyapunov construction. 
However, the existence of two purely imaginary 
eigenvalues means that the target state is no longer a 
hyperbolic fixed point but there is a center manifold 
of dimension two. Center manifold theory shows 
that the qualitative behavior near the fixed point is 
determined by the qualitative behavior of the flows 



on the center manifold [24]. Therefore, the next step 
is to determine the center manifold. For dimensions 
> 2 this is usually a non-trivial problem. However, 
if we can find an invariant manifold those tangent 
space at s d equals the tangent space of the center 
manifold, then this manifold is the center manifold. 
In our case solutions in the invariant set form a 
manifold diffeomorphic to the Bloch sphere for a 
qubit system, with the natural embedding 

P22 



p 



A*3 ft 



33 



1 



p 



A*3 ^33. 



which maps the state pd (or s d ) of the qutrit to the 
point s* d with p' d = diag(u>i, w 3 )/(wx + w 3 ), and 
the two tangent vectors of the center manifold at 
Pd to the two tangent vectors of the Bloch sphere 
at s* d . Thus this manifold is the required center 
manifold at pd (or s d ). On the center manifold 
Pd is a center with the nearby solutions cycling 
around it. The Hartman-Grobman theorem in center 
manifold theory proved by Carr [24] shows that 
all solutions outside E converge exponentially to 
solutions on the center manifold belonging to Sd, 
while the solutions actually converging to Sd only 
constitute a set of measure zero. Therefore, almost 
all solutions near p d converge to solutions on the 
center manifold other than p d and p d becomes no 
longer asymptotically stable (see Fig. [jjc)). 

B. H not strongly regular 

Let Hi fully connected but H not strongly reg- 
ular, e.g., 





"0 
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, H 1 = 


1 





1 










2lo 




1 


1 






Analogously to the section above, we can show 
that for a given stationary generic pd, the LaSalle 
invariant set forms a center manifold with the target 
state as a center. Hence, almost all trajectories 
near pd converge to other solutions on the center 
manifold and Sd is not asymptotically stable. 

VII. Concluding discussion 

We have studied a control design for tracking 
natural trajectories of generic quantum states based 
on the Hilbert- Schmidt distance as a Lyapunov 
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function. The analysis shows that the method is 
effective for generic density operators if and only 
if the invariant set is regular, i.e., contains only 
trajectories corresponding to critical points of the 
Lyapunov function. Since the Lyapunov function 
has exactly n\ isolated critical points for generic 
states, regularity of the invariant set in this case im- 
mediately implies that the target state or trajectory 
is isolated and thus locally asymptotically stable, 
but a detailed analysis shows that we have almost 
global convergence in this case. Although the set of 
states that do not converge to the target state is larger 
than previously asserted [10], for ideal systems it 
is still only a small subset of the state space, and 
for stationary target states we can show that it has 
measure zero. When the LaSalle invariant set is 
not regular in the other hand, the method not only 
becomes ineffective, and the target state ceases to be 
locally asymptotically stable, but a center manifold 
emerges around the target state, which exponentially 
attracts all trajectories, preventing convergence to 
the target state. 

The results follow from several steps. Computa- 
tion of both the LaSalle invariant set and the set of 
critical points of the Lyapunov function shows that 
a necessary condition for regularity of the invariant 
set is that the system Hamiltonian satisfy certain 
rather strict conditions, effectively equivalent to 
controllability of the linearization. Further analysis 
shows that when we restrict our attention to generic 
states, the Lyapunov function is a Morse function 
with n\ isolated critical points and the target state 
as the unique global minimum in addition to a 
unique global maximum and n\ — 2 saddle points. 
The critical points of the Lyapunov function further 
correspond to fixed points or critical trajectories of 
the dynamical system. If the dynamical system were 
a gradient flow of the Lyapunov function this would 
allow us to almost immediately infer almost global 
convergence to the target states. As this is not the 
case we must analyze the linearization of the dy- 
namics about the critical points and show that they 
are hyperbolic. We do this rigorously for stationary 
target states, where the analysis shows that the n\ 
critical points of the Lyapunov function are indeed 
hyperbolic fixed points of the dynamical system if 
the system Hamiltonian ideal. For stationary target 
states this condition also implies for regularity of 
the invariant set, and as the Lyapunov function is 
a Morse function in our case, we can use it to 



compute the dimensions of the stable and unstable 
manifolds at each of the hyperbolic critical points 
of the dynamical system. This shows that all critical 
points except the target state have stable manifolds 
of dimensions less than the state space and allows 
us to conclude that almost all initial states will 
converge to the target state in this case. The flipside 
of this analysis is that the target state ceases to be a 
hyperbolic fixed point of the dynamical system if the 
system Hamiltonian is no longer ideal, and in this 
case a center manifold emerges around the target 
state, which exponentially attracts all trajectories. 
For non-stationary target states the method can fail 
even if the system Hamiltonian is ideal, for target 
states that give rise to a non-regular invariant set, but 
we also show that such target states are a measure- 
zero subset of the state space. 
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Appendix 
A. Lie algebra generators 

A standard basis for the Lie algebra su(n) is given 
by {\ke, Afc^, \k} for 1 < k < £ < n, where 

Afc = i(e k k - e fe +i,fc+i) (25a) 
Am = i(eu + hk) (25b) 
\i = (he - hk) (25c) 

and the (k,£) th entry of the matrix e mn is <5 fcm <5 fa , 
and % = ^/— T. We have the useful identities 

Tr(A^Afc/^) = Tr(Aj^Afc'^) = —25kk'Sw (26a) 
Tr(X ke X k 'e) = (26b) 
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and for any diagonal matrix D = Y%=i dk&kk 

[D, X k ) = 0, (27a) 
[D, Xkt] = +i(d k - di)X k e, (27b) 
[D, X ke ] = -i{d k - d e )X ke , (27c) 

The basis ( [25] ) is not orthonormal but we can 
define an equivalent orthonormal basis {a m }mZi 
for su(n) by normalizing the n 2 — n non-Cartan 
generators 4^X k t and -j^\ k i, and defining the 



V2' 



n - 



1 



orthonormal generators for the Cartan subalgebra 

re r+liT+1 ) for 



cr„2 



n^—n+r 



: , [r(r + 1)] -i/2 ( ELl e ss 
n — 1. 



B. Lemmas in the proof of Theorem V.2 

Lemma A.l. Df(s ) vanishes on the subspace Sc- 

Proof: To show that Df(s )s = for all s G 
Sc, it suffices to show that A s = and s^Ais = 
for s G Sc. s G Sc corresponds to density operators 
p G iC, i.e., p diagonal. As A s is the Bloch 
vector associated with [— iH , p], —iH is diagonal 
and diagonal matrices commute, [— iH , p] = and 
A s = follows immediately. To establish the 
second part, note that for ip G C and — iU\ G T, 
we have [— iHi,ip) G T, or [—iHi,p] G iT, and 
Ais G St- Since p d is diagonal, i.e., s, d G S c _L Sr. 
we have s^A^s = for s*G S c . ■ 

Lemma A.2. The restriction of -D/(s*o) to St is 
well-defined and its matrix representation B has 
n 2 — n non-zero eigenvalues. 

Proof: Since we already know that Sc is in the 
kernel of Df(so), it suffices to show that the image 
of Df(s ) is contained in St, i.e., Df(s )sE St- 

D f (s )s = A s + A^o s d T Ais 
= A s + (sd T Ais) Ais 

shows that it suffices to show that A s G St and 
AiSq G St- Both relations follow from the fact that 
the commutator of a Cartan element and a non- 
Cartan element of the Lie algebra su(n) is always in 
the non-Cartan algebra T, and thus [— iH , p] G iT 
since — iH G C, and [— iH\, pd] G iT since ipd G C. 
Therefore, the restriction B : Sq — > St of Df(so)s 
is well defined. 

Furthermore, the restriction of A to St is a 
block-diagonal matrix B = di&g(A^'^) with 

_ . . [ ^ 



The restriction w of Ais to SV is a column vector 

of length n{n — 1) con- 



,^1,2) ^1,3) 



U 



(n— l,fi)' s 



sisting of 77,(72 — l)/2 elementary parts 



(M) 



A. 



(fc)r(f) 



v/2 



(28) 



for fc = 1, . . . , n — 1 and £ = k + 1, . . . , n. Similarly, 
let v be the restriction of Ai$d to St- Then ?7 = 



(1,2). _ _ _ . -(n-l,n) 



with £A fc '^ as in Eq. (28) and r 



the identity permutation. 

Thus the restriction of Df(s ) to the subspace 
St is B = B — uv T . Since u k £ ^ for all k, £ by 
regularity of H , we have det(B ) = Y[k,i w ki 7^ 0, 
i.e., B invertible, and det(B) = det(B — uv r ) = 
(1 — if Bq 1 u) det(B ) by the matrix determinant 
lemma [23]. Bq 1 is block-diagonal with blocks 



(M)i-i 



1 







Thus t^-B, 



it 



E k A^ k " 



"0 


-1" 






1 


_ 







■\T(j(k,£)fi(k,e) vanishes as 



0, Vfc, 



Hence, det(£>) = det(-Bn) 7^ 0, and the restriction 
of Df(s ) to St has only non-zero eigenvalues. ■ 

Lemma A.3. If if3 is a purely imaginary eigenvalue 
of B then it must be an eigenvalue of Bq, i.e., 
iP = HuJki for some (k, £), and either the asso- 
ciated eigenvector e must be an eigenvector of B 
with the same eigenvalue, or the restriction of A\s$ 
to the (k, £) subspace must vanish. 

Proof: If 27 is not an eigenvalue of B then 
(£>o — i(3I) is invertible and by the matrix determi- 
nant lemma 

= det{B - uv T - i(3I) 
= det((B - if3I) - uv T ) 
= (1 - v T (B - ipi)- 1 ^ det(B - i(3I). 

Since det(I? — i(3I) we must therefore have 



v T {B -i(5iy 



u 



1. 



Noting {B — i/3I) 1 is block-diagonal with blocks 



-i/3 
u k t 



(3 2 



-u k i 
-ip 



-iP -LO ki 
Um -iP 

»(6 M 



(29) 



—i 
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for all k, I, this leads to 



v T (B 
-i/3 



E 

kl 



k,£ 



M 



% • 



Since all terms in the sum are real this is a contra- 
diction. Thus if i/3 is a purely imaginary eigenvalue 
of B then it must be an eigenvalue of B . 

Since the spectrum of B is {±.iujkt}, this means 
i(3 = rLiuj k i for some (k,£). Without loss of gener- 
ality assume 7 = Ui 2 > and let e = x + iy be the 
associated eigenvector of B. Then 

uv T )(x + iy) = itu 12 (x + iy), (30) 



Be = (B 
which is equivalent to 

(B 



uv )x 



->-<T\ - 

uv )y 



Multiplying (31:)) by —uj\ 2 B y 
d3lk) 



0J12S. 
1 





(31a) 
(31b) 

and adding it to 



Bqx — uv x + luuBq uv y = —ui 12 B x 



Eq. (29) shows that 



12 



[B, 



M- 1 



5, 



(1,2) 



, i.e., on 



the T12 subspace the underlined terms above cancel, 
and thus the first two rows of the above system of 
equations are 

"0 -11 
1 



"1 

u 2 



(v r x) 



u 2 



v y). 



-C 2 Ui 



If v x 7^ then the last equation gives Ui 
and u 2 = —c 2 u 2 with c = v T y/v T x, which can only 
be satisfied if u x =u 2 = 0. Similarly if vFy 7^ 0. If 
iFx = v^y = then we have Be = B e = iuj\ 2 e, 
implying that e is an eigenvector of B associated 
with iuj 12 . ■ 



C. Proof of Lemma V.4 



Lemma A.4. Given a generic Pd(t), the LaSalle 
invariant set E is irregular if and only if 
rank A(sd(0)) < n 2 — n. 

Proof: It suffices to show that if rank A(sd) = 
n 2 — n, then for any p such that [p, Pd(0)} diagonal, 
we have [p, pd(0)} = 0. If this is true then for any 
(Ph P2) € E with diagonal commutator, there exists 
some t such that p 2 = e lH ° to pd(0)e~ lH ° to and since 
[pi,p 2 ] is diagonal, [ e ~ iH ^ Pl e iHoto , p d (0)} is also 
diagonal, hence equal to zero and [pi, p 2 ) = 0. 



Let p 2 = pd(0). First we show that the kernel of 
A(s 2 ) has dimension n — 1 and thus rank A(s 2 ) < 



n 



n. In this case rank A(sd) 



'ri- 



ii 



rank implies that the remaining n— 1 rows of 
A(s 2 ) are linear combinations of the rows of A(s 2 ) 
and therefore A(s 2 )s\ = implies A(s 2 )si = 0, or 
\PUP2] = 0. 

In order to show that the kernel of A (s 2 ) 
has dimension n — 1, we recall that if p 2 = 
U diag(wi, . . . , WnjW for some U G §U(n) then 
[Pi, P2) = for all pi = U diag(u; r( i), . . . , w T{n) )U\ 
where r is a permutation of {l,...,n}. If the 
Wk > are distinct then these pi's span at least 
a subspace of dimension n since the determinant of 
the circulant matrix 



C 



W 1 W 2 ... W n _i w n 
W 2 W 3 ... W n Wi 



W n _i W n ... W n _ 3 W n _ 2 
. W n Wi ... W n - 2 W n -!. 

is non-zero, and hence its columns are linearly 
independent and span the n-dimensional subspace 
of diagonal matrices. If the Wk are distinct then the 
kernel cannot have dimension greater than n — 1 
since the pi can only span a subspace isomorphic 
to the set of diagonal matrices. Thus, the kernel 
of A (s 2 ) has dimension n — 1. (The dimension is 
reduced by one since we drop the projection of 
p onto the identity in the Bloch representation.) 
Similarly, we can prove if rank A ( Sd(0)) < n 2 — n, 
then E contains points with nonzero commutator. 
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